load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/ut_string.ncl"

begin
    a = addfile("/public/users/wrf/ncepdata/Reanalysis/2011/hgt.2010.nc", "r")
;    printVarSummary(a)

    ; read time
    time = a->time
    ntime = dimsizes(time)
    ; convert UTC time to CST time
    time = time + 8
    timeStr = ut_string(time, "%Y-%N-%D_%H")
    timeStr_output = ut_string(time, "%Y-%N-%D_%H:%M:%S")

    ; read hgt
    hgt_raw = a->hgt
    lon = a->lon
    lat = a->lat
    hgt = (/hgt_raw*hgt_raw@scale_factor+hgt_raw@add_offset/) ; Unit: m/s
    copy_VarAtts(hgt_raw, hgt)
    copy_VarCoords(hgt_raw, hgt)
    delete(hgt@scale_factor)
    delete(hgt@add_offset)
    printVarSummary(hgt)

    ; plotting
;    wks_type = "x11"
    wks_type = "png"
    wks_type@wkHeight = 900
    wks_type@wkWidth = 900

    ; Create Plot Resource
    res                     = True
    res@cnFillOn            = True
    res@gsnSpreadColors     = True
    res@gsnAddCyclic        = False
    res@mpDataBaseVersion       = "MediumRes"
    res@mpDataSetName           = "Earth..4"
    res@mpGeophysicalLineColor      = "Black"
    res@mpNationalLineColor         = "Black"
    res@mpUSStateLineColor          = "Black"
    res@mpGridLineColor             = "Black"
    res@mpLimbLineColor             = "Black"
    res@mpPerimLineColor            = "Black"
    res@mpGeophysicalLineThicknessF = 2.0
    res@mpGridLineThicknessF        = 1.0
    res@mpLimbLineThicknessF        = 1.0
    res@mpNationalLineThicknessF    = 2.0
    res@mpUSStateLineThicknessF     = 1.0
    res@mpOutlineBoundarySets       = (/"National","Continental", "Geophysical"/)   ; turn on country boundaries
    res@mpOutlineSpecifiers = "China:states"

    ; crop China Region
    minLat = 15
    maxLat = 50
    minLon = 70
    maxLon = 130
    hgt_China = hgt(:,{500}, {minLat:maxLat}, {minLon:maxLon})
    printVarSummary(hgt_China)
    res@mpMinLatF           = minLat
    res@mpMaxLatF           = maxLat
    res@mpMinLonF           = minLon
    res@mpMaxLonF           = maxLon
    res@cnLevelSelectionMode = "ManualLevels"
    res@cnLevelSpacingF     = 60.
    res@cnMaxLevelValF      = 5900.
    res@cnMinLevelValF      = 5200.


    ; for Gansu North Grid
    hgt_Gansu_North = hgt_China(:, {40.7140}, {96.2402})
    printVarSummary(hgt_Gansu_North)

    ; for Gansu South Grid
    hgt_Gansu_South = hgt_China(:, {34.3525}, {104.8754})
    printVarSummary(hgt_Gansu_South)

    ; output Gansu Pressure to TXT file
    line_array = new(ntime + 1, "string")
    l_cnt = 0

    line_array(l_cnt) = "Time, Gansu North(40, 95), Gansu South(35, 105)"
    l_cnt = l_cnt + 1
    do it = 0, ntime-1, 1
        line_array(l_cnt) = timeStr_output(it) + ", " + \
            sprintf("%10.4f", hgt_Gansu_North(it)) + ", " + \
            sprintf("%10.4f", hgt_Gansu_South(it))
        l_cnt = l_cnt + 1
    end do
    hgt_txt_filename = "~/Data/NCEP/hgt_500hPa_2010.txt"
    asciiwrite(hgt_txt_filename, line_array)

    do it = 0, ntime-1, 1
        output_path = "~/Data/NCEP/figures/500hPa_HGT/"
        system("mkdir -p " + output_path)
        output_filename = output_path + "NCEP_500hPa_HGT_" + timeStr(it)
        print(output_filename)
        wks = gsn_open_wks(wks_type, output_filename)
        gsn_define_colormap(wks, "BlAqGrYeOrRe")

        hgt_China_it = hgt_China(it, :, :)
        hgt_China_it@long_name = "500hPa Height at " + timeStr(it)
;        printVarSummary(pres_China_it)
;        printVarSummary(pres_China)

        plot = gsn_csm_contour_map(wks, hgt_China_it, res)
        delete(wks)
    end do

end


